
#include "gsl_ext_function.h"
#include "gsl_ext_multimin.h"
#include "gsl_ext_testfunctions.h"
#include "GSLFunctionWrapper.h"
#include "GSLMinimizerWrapper.h"
#include "StoppingCriterion.h"
#include "SymbolicExpressionGenerator.h"
#include "OctaveArgParser.h"
#include "OctaveUtils.h"

#include <stdexcept>

GSLFunctionWrapper::Type OctaveArgParser::getGSLFunctionType(const string &minimizerName)
{
  if(minimizerName == "gsl_nmsimplex" || 
     minimizerName == "gsl_ext_lrwwsimplex")
    return GSLFunctionWrapper::F;
  else if(minimizerName == "gsl_ext_mnewton")
    return GSLFunctionWrapper::FD2F;
  else
    return GSLFunctionWrapper::FDF;
}

GSLMinimizerWrapper *OctaveArgParser::getMinimizerInstance(const string &name,
                                                           const Octave_map &params,
                                                           const Octave_map &linminParams,
                                                           int n,
                                                           double &stepSize,
                                                           double &tol)
{
  stepSize = 1.0;
  tol = 0.1;

  if(linminParams.contains("stepsize"))
    stepSize = linminParams.seek("stepsize")->second(0).double_value();

  if(linminParams.contains("tol"))
    tol = linminParams.seek("tol")->second(0).double_value();

  const gsl_multimin_fminimizer_type        *T_f    = NULL;
  const gsl_multimin_fdfminimizer_type      *T_fdf  = NULL;
  const gsl_ext_multimin_fd2fminimizer_type *T_fd2f = NULL;

  if(name == "gsl_ext_bfgs_hess_f")
    T_fdf = gsl_ext_multimin_fdfminimizer_bfgs_hess_f;
  else if(name == "gsl_ext_bfgs_f")
    T_fdf = gsl_ext_multimin_fdfminimizer_bfgs_f;
  else if(name == "gsl_ext_bfgs_hess_mt")
    T_fdf = gsl_ext_multimin_fdfminimizer_bfgs_hess_mt;
  else if(name == "gsl_ext_bfgs_mt")
    T_fdf = gsl_ext_multimin_fdfminimizer_bfgs_mt;
  else if(name == "gsl_ext_conjgrad_fr_mt")
    T_fdf = gsl_ext_multimin_fdfminimizer_conjgrad_fr_mt;
  else if(name == "gsl_ext_conjgrad_pr_mt")
    T_fdf = gsl_ext_multimin_fdfminimizer_conjgrad_pr_mt;
  else if(name == "gsl_ext_lrwwsimplex")
    T_f = gsl_ext_multimin_fminimizer_lrwwsimplex;
  else if(name == "gsl_ext_mnewton")
    T_fd2f = gsl_ext_multimin_fd2fminimizer_mnewton;
  else if(name == "gsl_conjugate_fr")
    T_fdf = gsl_multimin_fdfminimizer_conjugate_fr;
  else if(name == "gsl_conjugate_pr")
    T_fdf = gsl_multimin_fdfminimizer_conjugate_pr;
  else if(name == "gsl_nmsimplex")
    T_f = gsl_multimin_fminimizer_nmsimplex;
  else if(name == "gsl_steepest_descent")
    T_fdf = gsl_multimin_fdfminimizer_steepest_descent;
  else if(name == "gsl_vector_bfgs")
    T_fdf = gsl_multimin_fdfminimizer_vector_bfgs;
  else if(name == "gsl_vector_bfgs2")
    T_fdf = gsl_multimin_fdfminimizer_vector_bfgs2;
  else
    throw invalid_argument("invalid minimization algorithm name");

  if(T_f != NULL)
  {
    gsl_multimin_fminimizer *s_gsl = 
      gsl_multimin_fminimizer_alloc(T_f, n);
    if(s_gsl == NULL)
      throw runtime_error("failed to allocate memory");
    return new GSLMinimizerWrapper_f(s_gsl);
  }
  else if(T_fdf != NULL)
  {
    gsl_multimin_fdfminimizer *s_gsl = 
      gsl_multimin_fdfminimizer_alloc(T_fdf, n);
    if(s_gsl == NULL)
      throw runtime_error("failed to allocate memory");
    return new GSLMinimizerWrapper_fdf(s_gsl);
  }
  else if(T_fd2f != NULL)
  {
    gsl_ext_multimin_fd2fminimizer *s_gsl = 
      gsl_ext_multimin_fd2fminimizer_alloc(T_fd2f, n);
    if(s_gsl == NULL)
      throw runtime_error("failed to allocate memory");
    return new GSLMinimizerWrapper_fd2f(s_gsl);
  }
  else
    return NULL;
}

GSLFunctionWrapper *OctaveArgParser::getFunctionInstance(const string &name,
                                                         const Octave_map &params,
                                                         const Octave_map &evalParams,
                                                         GSLFunctionWrapper::Type type,
                                                         gsl_ext_testfunction_params *tfParams)
{
  string fEvalTypeStr, gEvalTypeStr, HEvalTypeStr;
  gsl_ext_func_eval_type fEvalType;
  gsl_ext_deriv_eval_type gEvalType;
  gsl_ext_deriv_eval_type HEvalType;
  tfParams->m = tfParams->n = 0;
  string expression;
	
  if(params.contains("m"))
    tfParams->m = params.seek("m")->second(0).int_value();
  
  if(params.contains("n"))
    tfParams->n = params.seek("n")->second(0).int_value();
	
  if(evalParams.contains("f"))
    fEvalTypeStr = evalParams.seek("f")->second(0).string_value();
	
  if(evalParams.contains("g"))
    gEvalTypeStr = evalParams.seek("g")->second(0).string_value();
	
  if(evalParams.contains("H"))
    HEvalTypeStr = evalParams.seek("H")->second(0).string_value();
	
  if(fEvalTypeStr == "sym")
  {
    fEvalType = GSL_EXT_FUN_EVAL_SYM;
    
    if(gEvalTypeStr == "sym")
      gEvalType = GSL_EXT_DERIV_EVAL_SYM;
    else if(gEvalTypeStr == "fd")
      gEvalType = GSL_EXT_DERIV_EVAL_FD;
    
    if(HEvalTypeStr == "sym")
      HEvalType = GSL_EXT_DERIV_EVAL_SYM;
    else if(HEvalTypeStr == "fd")
      HEvalType = GSL_EXT_DERIV_EVAL_FD;
    
    expression = 
      SymbolicExpressionGenerator::generateExpression(name,
                                                      tfParams->m,
                                                      tfParams->n);
    
    if(expression != "")
    {
      GSLFunctionWrapper *f;
      if(type == GSLFunctionWrapper::F)
      {
        gsl_multimin_function *f_gsl = 
          gsl_ext_multimin_function_alloc((char *)expression.c_str(), 
                                        fEvalType, 
                                        tfParams);
        if(f_gsl == NULL)
          throw invalid_argument("invalid objective function or parameters");
        f = GSLFunctionWrapperFactory::getInstance(f_gsl);
        return f;
      }
      else if(type == GSLFunctionWrapper::FDF)
      {
        gsl_multimin_function_fdf *f_gsl = 
          gsl_ext_multimin_function_fdf_alloc((char *)expression.c_str(), 
                                            fEvalType, 
                                            gEvalType, 
                                            tfParams);
        if(f_gsl == NULL)
          throw invalid_argument("invalid objective function or parameters");
        f = GSLFunctionWrapperFactory::getInstance(f_gsl);
        return f;
      }
      else
      {
        gsl_ext_multimin_function_fd2f *f_gsl = 
          gsl_ext_multimin_function_fd2f_alloc((char *)expression.c_str(), 
                                             fEvalType, 
                                             gEvalType, 
                                             HEvalType, 
                                             tfParams);
        if(f_gsl == NULL)
          throw invalid_argument("invalid objective function or parameters");
        f = GSLFunctionWrapperFactory::getInstance(f_gsl);
        return f;
      }

      if(f == NULL)
        throw invalid_argument("invalid objective function or parameters");
    }
  }
  else if(fEvalTypeStr == "precomp")
  {
    if(gEvalTypeStr != "" && gEvalTypeStr != "fd")
      throw invalid_argument("invalid evaluator combination");
    fEvalType = GSL_EXT_FUN_EVAL_PRECOMP;
    gEvalType = GSL_EXT_DERIV_EVAL_FD;
    HEvalType = GSL_EXT_DERIV_EVAL_FD;
  }
  else
    throw invalid_argument("invalid function evaluator type");
  
  if(type == GSLFunctionWrapper::F)
  {
    gsl_multimin_function *f = 
      gsl_ext_multimin_function_alloc((char *)name.c_str(), 
                                    fEvalType, 
                                    tfParams);
    if(f == NULL)
      throw invalid_argument("invalid objective function or parameters");
    
    return GSLFunctionWrapperFactory::getInstance(f);
  }
  else if(type == GSLFunctionWrapper::FDF)
  {
    gsl_multimin_function_fdf *f = 
      gsl_ext_multimin_function_fdf_alloc((char *)name.c_str(), 
                                        fEvalType, 
                                        gEvalType, 
                                        tfParams);
    if(f == NULL)
      throw invalid_argument("invalid objective function or parameters");
    
    return GSLFunctionWrapperFactory::getInstance(f);
  }
  else
  {
    gsl_ext_multimin_function_fd2f *f = 
      gsl_ext_multimin_function_fd2f_alloc((char *)name.c_str(), 
                                         fEvalType, 
                                         gEvalType, 
                                         HEvalType, 
                                         tfParams);
    if(f == NULL)
      throw invalid_argument("invalid objective function or parameters");
    
    return GSLFunctionWrapperFactory::getInstance(f);
  }
}

StoppingCriterion *OctaveArgParser::getStoppingCriterionInstance(const string &name,
                                                                 const Octave_map &params)
{
  double eps = 0.0;
  double fMin;
  int n = 0;
  gsl_vector *xMin = NULL;
  
  if(name != "maxnumiter" && params.contains("eps"))
  {
    octave_value eps_ = params.seek("eps")->second(0);
    if(eps_.is_real_scalar())
      eps = eps_.double_value();
    else
      throw invalid_argument("invalid threshold value type");
  }
  else if(name == "maxnumiter" && params.contains("n"))
  {
    octave_value n_ = params.seek("n")->second(0);
    if(n_.is_real_scalar())
      n = n_.int_value();
    else
      throw invalid_argument("invalid threshold value type");
  }
  else
    throw invalid_argument("missing threshold value for stopping criterion");
  
  if(name == "fdisttomin" || name == "fdisttominrel")
  {
    if(params.contains("fmin"))
    {
      octave_value fMin_ = params.seek("fmin")->second(0);
      if(fMin_.is_real_scalar())
        fMin = fMin_.double_value();
      else
        throw invalid_argument("invalid minimum function value");
    }
    else
      throw invalid_argument("missing minimum function value");
    
    if(name == "fdisttomin")
      return new FDistToMinTest(fMin, eps, false);
    else if(name == "fdisttominrel")
      return new FDistToMinTest(fMin, eps, true);
  }
  else if(name == "gradient")
    return new GradientNormTest(eps);
  else if(name == "maxnumiter")
    return new MaxNumIterTest(n);
  else if(name == "simplexsize")
    return new SimplexSizeTest(eps);
  else if(name == "xdiff")
    return new XDiffTest(eps);
  else if(name == "xdisttomin" || name == "xdisttominrel")
  {
    if(params.contains("xmin"))
    {
      octave_value xMin_ = params.seek("xmin")->second(0);
      if(xMin_.is_real_matrix())
        xMin = OctaveUtils::OctaveVectorToGSLVector(ColumnVector(xMin_.vector_value()));
      else
        throw invalid_argument("invalid function minimizer");
    }
    else
      throw invalid_argument("missing minimizer");
    
    StoppingCriterion *sc;
    
    if(name == "xdisttomin")
      sc = new XDistToMinTest(xMin, eps, false);
    else
      sc = new XDistToMinTest(xMin, eps, true);
    
    gsl_vector_free(xMin);
    return sc;
  }
  else
    throw invalid_argument("invalid stopping criterion name");
  
  return NULL;
}
